An immuno-epidemiological model with waning immunity after infection or vaccination

In epidemics, waning immunity is common after infection or vaccination of individuals. Immunity levels are highly heterogeneous and dynamic. This work presents an immuno-epidemiological model that captures the fundamental dynamic features of immunity acquisition and wane after infection or vaccination and analyzes mathematically its dynamical properties. The model consists of a system of first order partial differential equations, involving nonlinear integral terms and different transfer velocities. Structurally, the equation may be interpreted as a Fokker-Planck equation for a piecewise deterministic process. However, unlike the usual models, our equation involves nonlocal effects, representing the infectivity of the whole environment. This, together with the presence of different transfer velocities, makes the proved existence of a solution novel and nontrivial. In addition, the asymptotic behavior of the model is analyzed based on the obtained qualitative properties of the solution. An optimal control problem with objective function including the total number of deaths and costs of vaccination is explored. Numerical results describe the dynamic relationship between contact rates and optimal solutions. The approach can contribute to the understanding of the dynamics of immune responses at population level and may guide public health policies.


Introduction
In many infectious diseases, immunity acquired from infection wanes over time.The same holds for immune responses elicited by vaccines.Typical examples are influenza and COVID-19 where immunity wanes within a few months (Rambhia and Rambhia 2019;Goldberg et al. 2022).The durability of natural immunity and immune responses triggered by vaccines are crucial for decision making and interventions in public health.Antibodies seem to be the protective mechanism for these infections but often more specific immune responses such as specific T cell groups are needed to build up immunity and maintain immune memory.Immunity waning is highly heterogeneous in the population between individuals and changes over time (Lavine et al. 2021).
Several mathematical models have been developed to assess effectiveness and the possibility of waning immunity after infection or vaccination (Montalbán et al. 2022;Iyaniwura et al. 2023;Pell et al. 2022;Gosh et al. 2022;Domenech de Celles et al. 2022;Veliov and Widder 2016).To a lesser extent, the models investigated the optimal timing of vaccine administration, accounting for the waning immunity between seasons for infectious diseases such as influenza (Costantino et al. 2019).A population with heterogeneous immunity is considered in Montalbán et al. (2022).However, individual immunity is modeled as constant over time.In addition, Montalbán et al. (2022) consider no change in immunity levels due to previous infection or vaccination and do not study decision (control) aspects.Iyaniwura et al. (2023) used a distributed delay equations framework to describe the dynamics of waning immunity in a population with vaccine or natural infection induced immunity at an endemic stage.They performed a bifurcation analysis showing that waning immunity from natural infection influences the bifurcation type more than vaccine associated waning immunity.Furthermore, they derived a control reproduction number and showed the interplay between the decrease in immunity rate and the transmission rate of the pathogen.Similar approaches were used by Pell et al. (2022) and Gosh et al. (2022).Domenech de Celles et al. (2022) showed in a simulation study how immunological heterogeneity plays a role in determining the durability of vaccine protection.A model with heterogeneous dynamic immunity where sub-populations were structured with respect to the host immunity was developed and analysed by Veliov and Widder (2016).In all these cases investigation of control aspects was either not present or played a rather limited role.
In Sect. 2 of this study we propose a model that examines the dynamics of an infectious disease, taking into account the waning immunity following natural infection or vaccination.It is designed with the following key considerations: (1) The individuals are heterogeneous with respect to their (dynamic) immunity level; (2) After infection, individual immunity increases progressively until recovery begins; (3) With the onset of recovery, immunity starts to decrease over time; (4) infectiousness, susceptibility, mortality, and recovery rates depend on the individual immunity; (5) The infectiousness of the environment is represented by the aggregated infectiousness of infected individuals, weighted by activity level, as a share of the overall activity level of the population.
The model is formulated in terms of a system of Partial Differential Equations (PDEs) where the latter feature introduces a nonlocal effect in the form of a nonlinear term that incorporates integrals of state variables across the whole range of immunity levels.
From a mathematical perspective, the proposed model is challenging for the following reasons: (i) It consists of a system of first order PDEs (each of which is of size-structured type, see, e.g., Martcheva and Pilyugin (2006)) with different velocity fields, hence, with different characteristic lines.This creates a substantial problem in the analysis of the system, because a reformulation of the PDE system as an Ordinary Differential Equation (ODE) system in a closed form cannot be obtained.(ii) Because of the form of the non-local term describing the infectiousness of the environment, the Lipschitz constant of the equations may tend to infinity along the solution, which substantially complicates the proof of global existence.
As previously mentioned, several authors have investigated disease dynamics, considering factors such as waning immunity and the acquisition of immunity during infection or post-vaccination.Similar to our model is the work of White and Medley (1998), which involves equations with different transfer velocities.However, the authors focus on the formal steady-state equations, without examining the overall PDE system.Other studies, such as Rouderfer and Becker (1994), Barbarosa and Röst (2015), Ehrhardt et al. (2019), also consider first order PDEs, but either the velocity fiends are identical or a single PDE (together with ODEs) is involved.
Mathematical features of the model, such as the existence of a solution and the asymptotic behavior, are examined in Sect.3. The model is then extended in Sect. 4 to encompass effects of vaccination.Additionally, an optimal vaccination problem is formulated in Sect.4.2, which could potentially be utilized to design vaccine administration strategies.
Finally, in Sect. 5 numerical results are presented for several scenarios, which include the behavior of the epidemic with and without vaccination, as well as optimal vaccination policies.While mathematical properties of the PDE-system are analyzed in some detail, the optimal vaccination problem is analyzed only numerically.Here, analysis focuses on significant qualitative observations regarding the optimal vaccination policy and the corresponding evolution of the epidemic under optimal vaccination.
In the appendix of this work we prove global existence of a solution, even for more general systems than our particular model requires.The proof is not straightforward and may be of independent mathematical interest (see Sect. 3.2 for more explanations).

The basic model with dynamic immunity
To model the dynamics of the immunity over time, we use a function ω : R → [0, 1], whose value, ω(t), at time t is interpreted as the immunity level of an individual at time t.The larger this number ω(t), the higher is the immunity of the individual, which implies lower susceptibility and lower infectiousness.From an empirical point of view, the individual immunity level may be quantified, e.g., proportional the amount of antibodies per ml blood.
Throughout the paper, we assume that after an individual is infected, its immunity level increases until the time of recovery (Yaugel-Novoa et al. 2022).Therefore, we describe the evolution of the immunity level after infection at time τ using the equation where ξ ∈ [0, 1] is the initial immunity level at the time of infection τ .The function g : [0, 1] → [0, +∞) is assumed to be differentiable and to satisfy g(0) > 0, and g(1) = 0.The description we use for the immune response to the infectious agent corresponds to the way the immune response may be embedded in a dynamical system that describes the within-host pathogen dynamics of an infection (see Schuh et al. (2023)).In this approach the mounting and decline of the immune responses during an infection process were explicitly captured with an equation describing the overall immune capacity of the individual against the pathogen including the acqusition of an immunity level.In our approach we incorporate these immune response dynamics in an epidemiological model to capture the overall dynamics of the immunity at the population level.1 In the long run, the immunity level decreases (Yaugel-Novoa et al. 2022).Therefore, beginning with recovery immunity wanes over time and we describe its decrease by the equation where ξ ∈ [0, 1] is the immunity level at the time of recovery τ .It is assumed that While in a population the individual immunity levels change as described above, at any point in time the individuals in a population may have differing immunity levels, depending on their individual history with the disease.Therefore, the immunity status of a whole population can be modeled as a frequency distribution over the possible values of the immunity level, ω ∈ [0, 1].Note that ω here denotes just one possible value of the function ω(•).
Based on these considerations, we denote by S(t, ω) ≥ 0 and I (t, ω) ≥ 0 the size of the susceptible, respectively infected, population of immunity level ω at time t.Thus the total population N at time t is In this paper, we assume that the susceptibility and infectiousness of an individual depend only on its immunity level.Immunological memory may have been acquired through a history of previous exposure to the relevant pathogen through infection or vaccination.
We denote by c > 0 the contact rate of susceptible individuals, while the contact rate of infected individuals is represented by c I ∈ (0, c].In principle, the contact parameters can be extended to depend on ω, because people who know that they are well protected by immunity may have more contacts.Moreover, dependence on time may be used for the description of seasonal or other time dependent behaviour of the individuals.However, in this paper, we assume for simplicity that c and c I are constant. We model the rate of new infections by the expression cD(t)σ (ω)S (t, ω).This expresses the fact that the probability of an infection for each susceptible individual is proportional to its contact rate, its susceptibility and the infectiousness of the environment D(t).For one individual D can be considered as a stochastic process, depending on the actions of the individual and on the (random) infectiousness of its contacts.As we are finally interested in a model at the population level and consider classes of population groups with immunity level ω instead of individuals, we estimate the infectiousness of the environment as an average infectiousness.Infectious individuals of all immunity levels contribute with their infectiousness and at their contact rate to overall infectiousness.Here, the model is based on the mean field idea: for any contact, infectiousness and number of contactees are replaced by a population average.
In order to estimate the probability σ (ω)D(t) that an individual with immunity level ω is infected at time t, one has also to take into account the contact-adjusted size of the total population at time t.Therefore, under the assumption of weighted random mixing, the infectiousness of the environment in which susceptible individuals contact infected individuals is represented as Finally, the mortality rate of infected individuals is denoted by μ(ω), and the recovery rate from infection is denoted by ρ(ω).Both parameters are nonnegative functions depending on the current immunity level.In the present study we do not model in detail the demographic effects of birth and death rates, which may be even time varying or age-dependent.Neglecting their long-run effects, we basically assume that the birth rate and "natural" death rate are equal, and the same rates are effective for all relevant compartments.The mortality rate μ(ω) then represents the excess mortality due to the epidemic.
Based on these assumptions and the related notations, it is now possible to describe the time dependent dynamics of the classes of susceptible and infected individuals for different immunity levels in terms of a system of PDEs for the population sizes S and I of susceptible and infected individuals with varying immunity level.
(2.5) with initial conditions (S 0 and I 0 are initial data) and the zero flux boundary conditions Due to the assumptions f (0) = g(1) = 0 and f (1) < 0, g(0) > 0, the initial conditions and the zero-flux condition are equivalent to (2.7) Moreover, due to the meaning of ω, and for consistency of the initial and boundary conditions it is natural to assume that S 0 (1) = I 0 (0) = 0.In infectious disease epidemiology, modeling by size-structured systems is a well established approach, see e.g.Rouderfer and Becker (1994); White and Medley (1998); Martcheva and Pilyugin (2006); Barbarosa and Röst (2015); Veliov and Widder (2016); Ehrhardt et al. (2019).Each of the equations (2.4) and (2.5) is a standard size-structured equation.The equation represents the evolution of the concentration of a substance moving according to a given velocity field in presence of in-or outflow (the term on the right-hand side).It can be derived by the same (conservation of mass) argument as the advection equation (see e.g.Britton N.F. 1986) for a compressible gas.In contrast with the physical models, in population dynamics such equations usually contain non-local terms (the function D(t) in our case).Moreover, we deal with a system of two (later three) equations with different transfer velocities, which substantially complicates the analysis.
There is an alternative view of a system represented by equations (2.4), (2.4) (and subsequently, (4.5)), also widely used in mathematical biology.We discuss it in the remaining part of this section, where we normalize the population size such that Consequently, the compartment sizes, S and I (and V , in the next section) can be interpreted as proportions of the total population belonging to the respective subpopulations.
At the individual level we basically consider a Markovian stochastic process with hybrid state space: the state of any individual at any given time is characterized by their immune status and the compartment to which they belong at the time.The immune status takes continuous values, while the compartments are finite in number.In particular, compartments are the susceptible, the infected, subsequently also the vaccinated, and the dead individuals.Randomness is introduced only during the discrete transitions between compartments.In the intervals between these transitions, the continuous state progresses according to compartment-specific ordinary differential equations (2.2) and (2.1).Transitions take place in accordance with a Poisson process, with the deceased status serving as an absorbing state for all values of immunity level, and are governed by a time-and state-dependent infinitesimal generator If we assume for a moment that D(t) is a given function, the process can be considered as a piecewise deterministic process2 : in this case we may interpret the proportions S(t, ω), I (t, ω) as probabilities of being in the susceptible or infected state with immunity ω at time t.Then, equations (2.4)-(2.5),augmented by where G(t, ω) denotes the "probability" of being dead at time t (and having died at immunity level ω) is the Fokker-Planck (or Kolmogorov forward) equation of the process, see e.g., Annunziato and Borzi (2018).
The standard class of piecewise deterministic processes, sometimes also called "correlated random walk", was introduced by Davis (1984) and has been used extensively in theoretical biology, see, e.g., Rudnicki and Tyran-Kamińska (2000).Our model extends the standard piecewise deterministic case because D(t) is defined by (2.3), involving integrals over ω in a nonlinear way.Basically, the jump rate from S to I depends on the distribution of the whole population over all relevant compartments and all values of immunity levels at any point in time.

Notion of solution
The following assumptions hold throughout the paper.
For the (dummy) real numbers t, ω, d, s, i, denote (in relation to (2.4)-(2.5)) where the argument t is included for further use.For shortness we introduce the notations for γ ∈ , and for γ ∈ f and γ ∈ g , respectively.Figure 1 illustrates these notations.
The so-called "characteristic lines" staring from f cover , more precisely, the mapping {γ ∈ f , t ∈ [τ f (γ ), T ]} (γ , t) → (t, ω f [γ ](t)) ∈ is bijective.A similar fact applies to the characteristic lines emanating from g .Then any pair of continuous functions (S, I ) : → R 2 uniquely determines two family of functions of t parameterized by γ : These facts explain the following definition.
Definition 3.1 The pair of continuous functions S, I : → R is called a solution of system (2.3)-(2.7)if the functions z S and z I defined by (3.3)-(3.4)are absolutely continuous in t and satisfy the equations together with (2.3) and with initial conditions Equivalently, for any γ = (t, ω) ∈ it holds that (3.7) (3.8) Remark 3.1 We mention that if S and I are differentiable and satisfy equations (2.3)-(2.7) in the classical sense, then the corresponding representations on the characteristic lines solve equations (3.3), (3.4).Moreover, the representation (3.7), (3.8) is valid.The latter fact is not straightforward, however it can be directly checked using the identities and a few more similar identities that appear when plugging the expressions of S and I in (3.7)-(3.8)into (2.4)-(2.5).

Existence of a "smooth" solution
In this subsection we present a theorem claiming existence of a solution of system (2.3)-(2.7)which is regular enough to enable the subsequent analysis.Although the proof is based on the Banach contraction mapping theorem, it is not straightforward due to two reasons: (i) the Lipschitz constants of F S and F I may tend to infinity with the time due to the expression (2.3) for D, which makes the existence on [0, ∞) problematic; (ii) due to the involvement of different transfer velocity fields f and g, the system (2.4)-(2.5)cannot be reduced to a closed form ODE system along the characteristics; (iii) proving the non-negativity of the solution is not straightforward at all.Therefore, we present in the Appendix a detailed proof of the existence theorem formulated below.In fact, we even prove a more general theorem assuming a few properties of the functions F S and F I in (3.7)-(3.8)and not necessarily the specific form of (3.1)-(3.2).
Including more equations with different characteristic curves (such as the system with vaccination in the next section) does not change the proof.Since we allow dependence of the functions F S and F I on time in the proof, the presence of a control function in the equations is also covered by the existence theorem as proved in the Appendix.
The differentiability property in the claim of the theorem is crucial: it not only enables integration of equations (2.4)-(2.5)with respect to ω but also permits to interchange the order of integration and differentiation.This, in turn, implies that the aggregated compartment sizes can be represented as in the following corollary.
Then the aggregated compartment sizes Ŝ(t), Î (t) are given by (3.9) (3.10) In addition, the total number of individuals decreases according to Proof Integrate equations (6.1)-(2.5)over ω ∈ [0, 1] and exchange the order of integrals and differentials on the left hand side.This is possible due to the properties of S and I in Theorem 3.1.Apply the zero flux conditions (specified after (2.6)), to obtain expressions for d Ŝ(t) dt and d Î (t) dt .These derivatives exist for all t except for a finite number of points on every bounded set [0, T ].Finally, integrate over t ∈ [0, 1] to obtain equations (3.9)-(3.10).Equation (3.11) is then obtained by adding up.

Descend of the epidemics and basic reproduction numbers
The goal of this subsection is to obtain conditions under which the number of infected individuals decreases and converges to a disease-free state.We proceed under the assumption that the conditions stipulated in Theorem 3.1 are satisfied and introduce the notation Using the estimation which can easily be derived from the definition (2.3), we obtain by differentiating equation (3.10) (3.13) In two steps we eliminate the dependence of the estimation on the densities of S and I , focusing on the worst case: (3.15) Define the numbers Obviously λ t ≥ λ for any t ≥ 0. Thus we obtain the following proposition.

Proposition 3.3 At any time t, if the current normalized densities of susceptible and infected individuals, S(t,
•)/ Ŝ(t) and I (t, •)/ Î (t), satisfy the inequality λ t > 0 then the number of infected individuals strictly decreases at this time.Moreover, The next corollary claims that the susceptible population does not asymptotically extinct.Furthermore, it provides an estimate of the maximum population reduction attributable to the disease.

Corollary 3.4 Assume that λ > 0 and denote b
Proof From equation (3.9) and (3.12) we have The last inequality follows from We mention that in absence of disease, i.e.I 0 (•) = 0, the density of the susceptible individuals converges to the Dirac delta function concentrated at zero, irrespectively of the starting distribution S(t, ω).To formally show this fact we consider the purely susceptible population starting from S 0 (ω) with 1 0 S 0 (ω) dω = 1.Then S(t, ω) and I (t, ω) = 0 solves (2.4)-(2.5),where S is the solution of with S(0, ω) = S 0 (ω) and S(t, 1) = 0. Solving the last equation along the characteristics, after some elementary calculation we show that for every ϕ ∈ L ∞ (0, 1) and there is t(ε) such that for any t ≥ t(ε) we have In particular, S(t, ω) = 0 for a.e.ω ∈ [ε, 1] for such t.Now, take an arbitrary ϕ ∈ W 1,∞ (the domain is [0, 1], therefore W 1,∞ is the set of all Lipschitz continuous functions with the usual norm).We can estimate (using the obvious relation For t ≥ t(ε) the last term can be estimated by We can define the Dirac δ-function as an element of the dual space (W 1,∞ ) * , and we can also view S(t, •) as such.Then Basic reproduction numbers.Next, we investigate the related issue of basic reproduction number.For this, we consider a "small" portion of infected individuals, I 0 (ω), ω ∈ [0, 1], inserted in the susceptible population S 0 (ω), ω ∈ [0, 1], with 1 0 (S 0 (ω) + I 0 (ω)) dω = 1.This initially infected population changes over time due to recovery, mortality and increased immunity level, according to the equation with side conditions I (0, ω) = I 0 (ω), I (t, 0) = 0.The solution, call it I 0 (t, ω) (in the same sense as defined in the previous subsections), can be presented as where

Thus, abbreviating κ
Below we need the linearization of the function D in (2.3) with respect to I = I 0 at I = 0 and S 0 , which is c I 1 0 i(ω )I 0 (t, ω ) dω .Following the terminology in Diekmann and Heesterbeek (2000), we represent the "next generation" of infected individuals resulting from the "small" group of initially infected population I 0 (ω) as the solution of the equation Integrating in ω (see Corollary 3.2 and its proof), we obtain for Î (t) := 1 0 I (t, ω) dω the expression Changing the variable ω = ω g [0, ξ](t) and substituting I 0 with y[ξ ] we obtain the expression Since by a standard argument ∂ ∂ξ ω g [0, ξ](t) = e t 0 g (ω g [0,ξ ](s) ds , we obtain that Integrating in t we obtain the expression The left-hand side represents the total amount of secondary infections directly caused by I 0 , while the last multiplier on the right is the total amount of initially infected.Thus we can define the basic reproduction number of the disease, under the assumption that exact information about the ω-density of the initially infected and the susceptible population: We mention that the above improper integral may diverge.Clearly, a natural sufficient condition for convergence is that μ(1) + ρ(1) > 1.
Changing the variable t with ω = ω g [0, ξ](t) and the variable s with η = ω g [0, ξ](s) we obtain the equivalent formula If exact information about the ω-density of the initially infected individuals is not available, a worst case scenario is included in the definition, as it is usual in the considerations of the basic reproduction number (see e.g.Von den Driessche and Watmough 2008).If only the density of S 0 is known considering the worst case of gives the expression where This leads to the definition of basic reproduction number of the disease, under the assumption that exact information about the ω-density of susceptible population is known: Finally, if no exact information about the initial density of the susceptible population is available, then the worst case scenario gives the basic reproduction number Notice that one can estimate R 0 [S 0 (•)] and R 0 using the obvious inequality In particular, this gives . (3.17) For comparison, we mention that in the case of data independent of ω, the estimation of basic reproduction number (3.17) reduces to c I iσ ρ+μ , which coincides with the usual expression for the basic reproduction in the SIR model.

Modeling and optimization of vaccination
In this section we introduce a control variable that represents the vaccination efforts and consider a class of optimization problems for the vaccination policy.

Modelling vaccination
Including vaccination requires an extension of the basic model 2.3-2.7.We assume that only susceptible individuals are vaccinated.It is necessary then to model the act of vaccination together with the immunological behavior of vaccinated individuals.
We denote by v(t, ω) the vaccination rate applied to susceptible individuals of immunity level ω at time t.This means, that v(t, ω)S(t, ω) individuals of immunity level ω become vaccinated at time t.
The effect of vaccination on immunity is not immediate.Like newly infected individuals, vaccinated individuals gain immunity over time, until their immunity level reaches a maximum, possibly depending on the immunity level before vaccination.After that, the immunity level slowly decreases in the same way as that of all susceptible individuals with the same immunity level (Goel et al. 2021).
Therefore, we augment the model by adding a compartment, representing newly vaccinated individuals acquiring immunity after vaccination.Similarly as for susceptible and infected individuals, the size of this compartment is counted separately for each immunity level over time, and is denoted by V (t, ω).The process of acquiring immunity from vaccination (in a relatively short period after vaccination) is modeled in a similar way to the increase of immunity during illness, namely by the equation where ξ ∈ [0, 1] is the initial immunity level at the time of vaccination τ .The function h(ω) : [0, 1] → [0, ∞), represents how fast immunity builds up at the current immunity level ω.Presumably, it is a decreasing function, with h(0) > 0 and h(1) = 0, that is, with similar properties as the function g.
When reaching their individual maximum immunity level, newly vaccinated individuals leave the vaccinated compartment and are counted as susceptible individuals with the attained new immunity level.This means that their decrease in immunity will change in the same way as for susceptible individuals of the same immune level.The transition process from vaccinated to susceptible occurs with a rate r (ω), so that 1/r (ω) is the average duration of increase in immunity depending on the current immunity level.
Since vaccinated individuals behave in their activities as susceptible, the infectiousness of the environment, D(t), takes the form 2) The overall model with vaccination takes the form with initial conditions and boundary conditions Remark 4.1 As long as the immunity level of people is not measured in reality, the dependence of v on ω is an idealization.The time from last vaccination or from the last infection could be considered (as practiced in reality) as a proxy for the individual immunity level.It is also possible to consider the vaccination rate depending on time only: v = v(t).In this case one can formally substitute v(t, ω) = v(t) in the equations, which means equal vaccination rate for all ω.
Remark 4.2 Since V (0, ω) = 0, the formulae for the basic reproduction numbers remain the same as in the no-vaccination case.

Optimal vaccination policies
In this subsection we use the model involving the vaccination rate v(t, ω) to formulate an optimal control problem, which reflects the desire of acting in an rational way.
Literature considers a number of reasonable objectives for public health interventions, in particular vaccination, involving the burden on hospitals, the number of sick individuals, the cost of policy measures, as well as direct or indirect economic losses due to the disease or due to the policy measures, see e.g.Bloom et al. (2020); Caulkins et al. (2021).We focus on the following three objectives (to be minimized) posed on a fixed time-horizon [0, T ].
(i) The most important objective is expressed by the total number of deaths; on average this number represents also the total number of sick people, hence it also reflects the economic cost of the epidemics in absence of additional restrictive measures such as partial lock-down (not employed anymore in most of the countries after 2022).
(ii) Due to various reasons (religious, disbelief, feeling of violation of freedom, etc.) a part of the society is not willing to vaccinate; this is the reason because of which in several European countries the vaccination level is rather low.The decision makers, i.e. governments have to take into account the social tension created by compulsory vaccination and the resulting "social disutility".Consideration of disutility directly resulting from policy measures is not typical in the literature, although it is a factor that often strongly influences the real decision maker (especially at a political level).We refer to Bloom et al. (2020), Section 4.2 (after (4)), where social disutility of policy measures is involved in the objective function.
(iii) The cost of vaccination, which is perhaps less significant than the first two especially in public health emergency situations.
The first objective is clearly contradictory to the other two.Therefore, in the spirit of Pareto's approach to multi-criteria optimization problems, we define the weighted aggregated objective to be minimized as Here, α ≥ 0 and β ≥ 0 are weighting parameters.The choice of the above quadratic specifications of the disutility function is somewhat arbitrary, although the quadratic form is not needed in the subsequent numerical simulations.This choice reflects the fact that the social disutility marginally increases with the vaccination rate.The optimization (minimization) is subjected to the constraints (4.2)-(4.7)and the control constraint v(t, ω) ≥ 0.
In this paper, we do not deeply investigate issues as existence of an optimal solution, necessary optimality conditions, convergence of numerical methods, etc.However, due to the linear-convex form of the objective functional and the linearity of the equations (4.3)-(4.5)with respect to v, one may expect that an optimal solution exists and the optimal control is Lipschitz continuous.Although this is far not enough to claim convergence of our numerical approach, the results of the numerical experiments (see the next section) and the pertaining sensitivity analysis support such a claim.
The numerical approach we employ is the so-called direct method in optimal control, which consists of direct discretization of the equations and the objective functional in time and space (for ω), as briefly described in Subsection 5.1.

Numerical experiments
In the following section we provide several purely illustrative numerical experiments for the evolution of the model dynamics with and without vaccination.Moreover, we also analyze the impact of optimal vaccination policies among sub-populations with differing immunity level.

Numerical approximation
In order to obtain numerical solution to (4.2)-(4.7)we use the so called upwind scheme which is of first order accuracy, see LeVeque 2002.Consider for example equation (4.3), which can be written in the following way: with the initial/boundary conditions (4.6) and (4.7).Here, Z (t, ω) := S(t, ω), I (t, ω), V (t, ω) and F S is the right-hand side of (4.3).
In order to describe the numerical scheme, we define the uniform mesh ω j , j = 1, ..., M + 1, in the ω-dimension with step size ω = ω i+1 − ω i .Similarly, in the t-dimension we define the mesh t i , i = 1, ..., N + 1, with step size t = t i+1 − t i .The upwind scheme is represented by the following implicit recurrent equation: for i = 1, ..., N and j = 1, ..., M. From the boundary condition (4.7) we have that S(t i , ω M+1 ) = 0, for every grid point t i .
The scheme has to take into consideration also the sign of the functions f , g and h.For the equations (4.4)-(4.6)we have to change f with g or h, the numerator on the right-hand side to I (t i , ω j ) − I (t i , ω j−1 ) or V (t i , ω j ) − V (t i , ω j−1 ), and account for the zero boundary condition (4.7).
A necessary condition for convergence of the upwind scheme is the Courant-Friedrichs-Lewy condition (CFL), see Courant (1967).In our case this condition takes the form u t ω ≤ C < 1, where u = max ω { f (ω), g(ω), h(ω)} and C is the Courant number for the problem.We tested the scheme with various step sizes that satisfy the CFL condition in order to obtain reliable numerical results.

Model parameters -the baseline scenario
We outline the parameters selected for a baseline case, which will be used in the subsequent numerical analysis.In addition, the base case values will be varied to study the solution sensitivity.It is important to note that all values are selected for illustrative purposes only and do not correspond to a specific disease: while the current work presents a general model and analyzes some of its properties, substantial empirical work remains to be done in order to apply it to real-world data.
The initial distribution of the compartment sizes is chosen consistently with the boundary conditions provided in (2.7), respectively (4.7).For the distribution of the susceptible population at time zero, we take a linear function S 0 (ω) = 1.9(1 − ω), so that 95% of the population is susceptible at time t = 0.A parabola is chosen for modeling the distribution of the initial infected population, I 0 (ω) = 0.3 ω(1 − ω), which gives initial infected population 5% of the total.There are no vaccinated individuals at the beginning: V 0 (ω) = 0.
The parameters and functions for modeling contact rates, infectiousness, recovery and mortality are summarized in Table 1.
Table 1 also shows the specific functions f , g and h, which are specified as affine, such that the sign of their slope determines gain or loss of immunity over time.
In order to propose reasonable functional specifications, the following factors were taken into account.Susceptibility, infectiousness and mortality, σ, ι, and μ, decrease when immunity is higher.The specific form of the recovery rates ρ and r ensures that the average duration of recovery is shorter when immunity is higher; for zero immunity ρ(0) = 1/40, thus for ω = 0 the average duration of recovery is 40 days; for ω = 0.5 the average duration is about 6 days, for ω = 1 it is negligible.

Numerical results without vaccination
We start with numerical simulations of the baseline model, excluding vaccination.Figure 2 gives an overview of the development of the epidemic over the simulation horizon with two graphs.The first wave of the epidemic plus the emergence of a second wave is shown in Fig. 2a.Here, the numbers of individuals in the compartments of susceptible and infected individuals are aggregated over all immunity levels ω. 1 0 ωS(t, ω) dω/ 1 0 S(t, ω) dω, etc.In addition, the average immunity level of newly infected and newly recovered individuals is illustrated.In particular, one can observe an increase of the average immunity level in all compartments after the number of infected individuals peaks (compare with the first plot 2b).In the susceptible compartment the average immunity increases due to the inflow of recovered individuals with high immunity level, until the number of infected individuals becomes low enough.The magnitude of the immunity increase for a given initial immunity level mainly depends on the choice of the function g.
Figure 3 presents the evolution of the densities of the newly infected individuals and the newly recovered individuals.For each time t (on the horizontal axis), the related vertical cut pictures the normalized densities σ (•)S(t, •)/ 1 0 σ (ω)S(t, ω) dω and ρ(•)I (t, •)/ 1 0 ρ(ω)I (t, ω) dω.The bright yellow spots correspond to higher densities and the dark blue color corresponds to a low density.The immunity of the newly recovered individuals is much higher than that of the newly infected, which is to be expected, given that the immunity level increases during the infection.In order to show the long term behavior of the model, we simulate the baseline scenario without mortality, μ(ω) = 0, on a longer period of time, T = 16000 days.As depicted in Fig. 4a, the system exhibits an oscillatory behaviour with declining amplitude, which suggests convergence towards an endemic equilibrium.Figure 4b shows the decrease of the group of infected individuals when the contact rate c I of infected individuals has the lower value 0.1.For this value the condition in Proposition 3.3 is satisfied and we observe convergence to an disease-free equilibrium.

Simulations with constant vaccination
Before considering the vaccination rate v(t, ω) as a control variable, we analyze in the current subsection the effect of a constant and ω-independent vaccination rate, v(t, ω) ≡ v, on the evolution of the epidemics in the baseline case.
In Fig. 5a one can see the number of newly vaccinated individuals, v 1 0 S(t, ω) dω. Figure 5b presents the number of infected individuals for various vaccination rates v. Observe that for v(t, ω) = 0.2 the number of infected individuals approaches zero after the first wave.
In the next experiment we consider the results of vaccination depending on the initial time of implementation of the vaccine.This is done for three constant vaccination rates: v = 0.15, 0, 20, 0.25.Figure 6 shows the saved lives per vaccination (compared with the case without vaccination) depending the on initial time.It is remarkable that the three curves are strictly convex, which means that a delay in the implementation of a vaccine causes more deaths in earlier stages of the epidemics than in later stages.One reason for that is the increase of the average immunity level in the course of epidemics due to infections (manifestation of the herd immunity in the present model) and the marginal increase of immunity due to vaccination is smaller.

Optimal vaccination
Now, we analyze numerically the control problem, formulated in Subsection 4.2.The control function v is considered as piece-wise constant.We use the direct method of optimal control and direct discretization of equations (4.3)-(4.5)by the scheme described in Sect.5.1.Furthermore, the integral in (4.8) is approximated by a second-order quadrature formula.For solving the resulting mathematical programming problem we utilize the SQP method available as a function of the MATLAB™ Optimization Toolbox.The parameters of the baseline scenario are shown in Table 1.
Figure 7 shows the effect of optimal vaccination in the baseline case.In addition to the overall numbers of susceptible and infected individuals, the yellow line in Fig. 7a indicates the total number of individuals in the vaccinated group.These are individuals who are in the process of acquiring immunity due to the vaccination.As expected, the overall level of infected individuals is greatly reduced compared to the baseline case without vaccination in Fig. 2a.Surprisingly, the reason for that is not that the average immunity of the non-infected population is standingly higher than that of the non-vaccinated population.In Fig. 7b we see the average immunity of all compartments.Comparing the blue line on the previous plot 2b with that on 7b we see that they are quite similar, with the difference that the average immunity of the Fig. 7 Evolution of epidemiological population with optimal vaccination applied Fig. 8 Optimal Strategy for the baseline scenario: Administered vaccines, i.e. newly vaccinated individuals and comparison with the average immunity level of the susceptible group vaccinated population increases to its maximum substantially earlier than that of the non-vaccinated population.Thus the immunity level of the infected population is much higher exactly in the expansion face of the epidemics which leads to less infections.Later on, the immunity level of the non-infected population catches up due to the higher herd immunity.
The optimal vaccination policy is analyzed in Fig. 8. Figure 8a shows the optimal total number of newly vaccinated individuals.It can be seen that the main effort should be concentrated immediately at the beginning in order to boost immunity in the population, while only a smaller effort is dedicated later-on in order to replace the waning immunity and to maintain a low level of infections.Due to the finite horizon of the optimal control, the vaccination is terminates before the end of the simulation period.
Because we allow vaccination efforts to depend on the immunity level, optimal vaccination policy is not only a matter of the overall level and timing of vaccination.Figure 8b shows the distribution of the application of vaccines to individuals with differing immunity level over time.Again, the abrupt change at the end of the horizon is due to the stop of vaccination.It can be seen that, at the beginning -when the overall effort is high, vaccination tends to be given to individuals with lower immunity level than the average immunity level in the susceptible group.The levels of immunity of vaccinated individuals then catch up with the average immunity level of susceptible individuals after around 150 days, and then follow the general decrease of the average level of immunity.
Dependence on the time-horizon and model predictive control.In practice, any vaccination policy has to be revised after some time to catch up with new information.In particular, the improved medical understanding of the disease, changes in the death rates, new variants of the pathogens may emerge and enhanced vaccines may be developed.In terms of control, a new optimization is done after some time with updated information, which is known as Model Predictive Control in the literature.The question arise, whether our model is suitable for such revisions.
In order to apply Model Predictive Control, we may solve the optimization problem on a relatively short horizon, e.g.t = 400 days in the baseline case, apply the obtained solution during an even shorter time horizon, say 70 days, then update the model parameters and the current real state of the epidemic, solve the problem on the next 400 days horizon, and so on.
Such an approach only works well if the results with different planning horizons do not vary too much over the shorter time period (here 70 days).This is tested in Fig. 9, which shows the dependence of the optimal vaccination policy on the chosen time horizon [0, t] on which the optimization problem is solved.The plot on the left indicates that the total number of optimally vaccinated people is practically independent of the time horizon over the first 70 to 100 days.More relevant, the same applies (even on a longer horizon) to the aggregated (in ω) vaccination effort (the right Fig. 9).So, it seems to be reasonable to apply Model Predictive Control.
Trade-off between the objectives.The objective (4.8) puts together three objectives: the number of deaths and two kinds of vaccination costs.All these goals are relevant for decision-making in public health, but also are contradictory.In order to formulate the overall objective, weights are applied, modeling the relative importance of the individual objectives.
Although the analysis of the optimization problem so far has been focused on the analysis of one ("standard") choice of these weights (see the baseline scenario defined in Table 1), it is possible to go deeper by analyzing the efficient frontier (or Pareto frontier).With conflicting goals, vaccination strategies can be compared by showing the vectors of their respective partial objectives (deaths, number of vaccinations, social cost) in a plot.Iterating over the possible weight combinations (α, β) and plotting the values of partial goals, one gets the Pareto frontier.From the viewpoint of a decision maker, the efficient frontier depicts those combinations of conflicting goals that are achievable at the best.
Having three objectives, the efficient frontier is a 2-dimensional manifold in the 3dimensional space.In order to show a two-dimensional picture, we vary only the weight of the administration cost, β, and plot the Pareto curve, holding fixed the parameter α as in the baseline case (the weight of the deaths is fixed to 1 by normalization of the overall objective function).Varying the weight β in the range from 0.001 to 0.02 and calculating the corresponding optimal vaccination policy, we obtain (a part of) the efficient Pareto curve.
We consider three cases for the vaccination control v: (i) the control may depend on time and immunity level as in the previous considerations, that is, v = v(t, ω); (ii) the control depends only on time, v = v(t), and the vaccination is uniformly distributed over individuals of different ω; (iii) the control is constant across time and immunity levels.The first case represents an idealized situation in which full information about the immunity level is required, while in the other two scenarios such information is not needed.This setup allows to quantify the effects of available information on the objective values, similar to the concept of the value of information in stochastic optimization.
Fig. 10a shows the efficient frontiers for the three scenarios, and also depicts the Summarizing, the blue curve in Fig. 10a represents the optimal vaccination administration cost (in term of the number of vaccinations) versus the optimal achievable total percentage of deaths for case (i).For any point (b, d) on this efficient frontier, d is the minimal percentage of deaths that can be achieved by vaccination budget b.
Vice versa, if it is intended to limit the number of deaths to d%, then at least b are the necessary vaccination costs.By the strict convexity of the Pareto curve, the less is the number of deaths, the more costly it becomes to obtain any additional reduction of the of deaths.Similar explanation applies to the other two control scenarios.Note that the baseline scenario with no vaccination at all is efficient for all three cases.It lies on the point, where all three curves touch the ordinate.Fig. 10b shows the overall control strategy (aggregated over ω in case (i)) for the three baseline cases, related to the red points in Fig. 10a.It suggests that, compared to the other two scenarios, more people should be vaccinated at the beginning of the epidemics, if information on the distribution of ω is available (control scenario (i)).The value of information (and of capacity to act) is demonstrated by the mutual positioning of the three curves in Fig. 10a: the Pareto frontiers of cases with more information or capacity to act, clearly dominate the other curves.

Comparative analysis
We stay with the baseline case and the optimal control problem, considered so far, and analyze the effect of parameter changes on the objective value.As contact rate plays a significant role in pathogen transmission, we vary the baseline values for c = 8 and c I = 3, by multiplying these parameters by a factor ranging from 0.8 to 1.25.For each parameter value we calculate the optimal vaccination policy and plot the corresponding optimal percent of deaths in Fig. 11a and the remaining part of the optimal objective value in (4.8) (representing the disutility from, and cost of vaccination) in Fig. 11b.In both plots, the x-axis shows the values of the contact rate c.
In the graph in 11a the lowest mortality cost results from the contact rate c = 8.8 and in the second plot 11b the highest vaccination cost occurs for c = 8.4.One can observe that a lower contact rate does not necessarily imply lower mortality.When the contact rate increases, the optimal vaccination efforts also increase, but only up to a point.The vaccination leads to lower mortality rates, but we can also observe from figure 11b, that the vaccination costs for highest values of the contact rate is reduced.This fact can be explained by the effect of herd immunity.With higher contact rate more people obtain immunity from infection and the effect of vaccination is relatively smaller.In this sense, the vaccination and the herd immunity appear as substitutes when the contact rates are sufficiently high.

Discussion
In this study, we introduce and analyze an epidemiological model that explicitly incorporates the impact of waning immunity following infection or vaccination.The model differs substantially from previous approaches in the epidemiological literature (White and Medley 1998;Rouderfer and Becker 1994;Barbarosa and Röst 2015;Ehrhardt et al. 2019).It consists of a system of two PDEs of first order.When vaccination is considered, a third equation is added.The complexity of the model lies in its mathematical intricacies, primarily the inclusion of a nonlocal term, which necessitates the integration of state variables in a nonlinear fashion.This complexity is further amplified by the fact that the velocity fields are different in different equations.
A qualitative study of the model is provided that includes existence of a global solution, conditions for decay of the epidemics from a given state are obtained, and basic reproduction numbers under various information patterns.
In a further step, we introduce vaccination strategies and formulate an optimal control problem with three objectives: the total number of deaths, the social discomfort created by the pressure that people experience when the vaccination effort is high, and the direct costs of vaccination.Using plausible scenarios of vaccination, numerical results provide insights into the dynamics of the epidemiological populations involved, including the waning immunity with and without vaccination.With respect to the optimal vaccination strategy the model provides insights into the influence of different factors on the optimal policy and performance.An interesting fact, for example, is that vaccination efforts and herd immunity act in a certain sense as substitutes: above a threshold value of the contact rate, further increase of the contact rate leads to lower vaccination rate.Below this threshold value, increase of the contact rate leads to increase of vaccination rate.In addition, we determine the efficient frontier between vaccine administration costs (direct and indirect) and the number of deaths, using three different control settings: optimal control policy which is independent of time and immune level, optimal policy that depends only on the time, and optimal policy depending on time and immunity level.
Although the model has some striking features such as the description and coupling of the elicitation of immune responses to the epidemiological process it can only be considered as a first step towards a more detailed description of immune responses and their waning over time.In particular, the model covers the immune response in a broad sense, without differentiating between innate and adaptive responses and their antibody and cellular branches.It does not account for their unique dynamical attributes, such as the time lag between the two, and their varying intensities contingent on the infectious disease under investigation.There is still substantial clinical and epidemiological empirical work to be done to this type of model with relevant immuno-epidemiological data.Despite its limitations, the model demonstrates that a mathematical representation of these dynamic processes is feasible.This could pave the way for a deeper comprehension of the processes in question and the assessment of related interventions.

Appendix: Proof of Theorem 3.1
Since the horizon T may change in the subsequent considerations, at some places we use the notation T := [0, T ] × [0, 1].The space of all continuous functions from a set X ⊂ R n to R is denoted by C(X ), with the usual norm denoted by . The spaces L 1 (0, T ) and L ∞ (0, T ) are defined as usual.Further, L(X ) ⊂ C(X ) is the subspace of all Lipschitz continuous functions with Lip(x) denoting the (minimal) Lipschitz constant of x ∈ L(X ).We abbreviate The existence theorem presented below will be formulated in the terms of a general function F in the equations (3.7), (3.8), and with a general relation between the functions D and (S, I ).Namely, instead of equation (2.3) we set where D has the following properties: there exists constants a > 1 and Keeping in mind the specific form of the functions F S and F I in (3.1), (3.2), we assume in the general case that there exist constants M and L such that Moreover, the following property is fulfilled: for any d ≥ 0, (t, ω) ∈ F S (t, ω, d, 0, i) ≥ 0 ∀i ≥ 0, F I (t, ω, d, s, 0) ≥ 0 ∀s ≥ 0. (6.6) Theorem 6.1 Let the functions f , g, S 0 , I 0 satisfy the Standing Assumptions (at the beginning of Subsection 3.1).Let, in addition, the conditions (6.2)-(6.6)be fulfilled.
Then there exists T > 0, independent of the particular initial data (S 0 , I 0 ), such that the system (3.7),(3.8), (6.1) has a unique nonnegative Lipschitz continuous solution Z = (S, I ) on T , satisfying the inequality D(Z ) ≤ a − 1.
Let us fix the number T > 0 such that Notice that T does not depend on the initial data (S 0 , I 0 ).
and define the set (2.6) and (2.7) are satisfied . (6.9) On K T ,a we define the mapping F [D] as (6.10) (6.11)where γ = (t, ω) ∈ T .In the next three parts of the proof we shall prove that K T ,a is a nonempty complete metric space, F [D] maps K T ,a to K T ,a , and F [D] is contractive.
Due to the uniform Lipschitz property in the definition of the set K T ,a , it is a complete metric space in the metric induced by the norm in C( T ).

(Invariance of K
satisfies the side conditions in (2.6) and (2.7).The same applies to F I [D] (Z ).Fix an arbitrary Z ∈ K T ,a .Using (6.4), we have for any Then Thus F [D] (Z ) fulfills the growth condition in the definition of K T ,a .For any Then we split the above integrals into three parts (in each of the integrals only two parts may be non-degenerate).The integration in the first part is on an interval of length In view of (6.4) and the growth condition in the definition of K T ,a , the integrands can be mojorated by Ma(2e 2Ma Z 0 C(0,1) ).Then the sum of these integrals is at most 4Ma where we make use of (6.4), the growth condition and the Lipschitz property in the definition of K T ,a .Combining the obtained estimations and using the same estimations for F I [D] (Z ), we obtain that where in the last inequality we have used (6.7) and (6.8).This completes the proof of the invariance of K T ,a .

(Contractivity of F [D]
.) For Z , Z 1 ∈ K T ,a we have, using (6.5) and (6.7), According to the Banach contraction mapping theorem, for any function
For a fixed D as above, we shorten the notation Z where C := a max{L, M}.Combining the two cases we obtain that By the same argument, a similar inequality is fulfilled by q.Summing the two inequalities we obtain that Since p and q are Lipschitz continuous non-positive functions, we conclude that p(t) + q(t) = 0 on [t, T ], hence also on [0, T ].Then (6.14) implies that p(t) = 0 and similarly q(t) = 0.This proves the nonnegativity of S and I .
C , according to the growth condition in the definition of K T ,a .A similar inequality holds for I .Summing the two, and taking the supremum in ω ∈ [0, 1] on the left-hand side, we obtain that Using the Grünwal inequality we obtain that where L Z := 4Lae 2aL+2a M Z 0 C .This inequality gives the meaning of the Lipschitz property of Z [D].
6. (Proof of the existence claim in Theorem 6.1.)Define the set Then the solution Z [D] of (3.7)-(3.8),defined in point 4 of the proof, exists for every D ∈ N T and we may define Due to the properties of D and N T the latter is invariant with respect G. Apparently, it is a complete metric space.We shall show that the mapping G is contractive with respect to the norm D Thus G is a contraction on N T hence it has a fix point D * .Obviously the pair (Z * := Z [D * ], D * ) satisfies the system (3.7),(3.8), (6.1).Moreover, Z * ∈ K T ,a is a Lipschitz function, which completes the proof of the existence part of Theorem 6.1.

7.
(Proof of the uniqueness claim in Theorem 6.1.)Let for some T * > 0 system (3.7),(3.8), (6.1) has two nonnegative Lipschitz continuous solutions Z i = (S i , I i ), i = 1, 2, on T * , satisfying D(Z i ) ≤ a − 1.Without any restriction we may assume that Z 1 and Z 2 differ from each other for some arbitrary small t > 0 (otherwise we may compare these solutions starting at a later time chosen so that the solutions immediately decline from each other).Consider any T ∈ (0, T * ] (to be fixed later).Denote D i := D(Z i ).Then |D i (t)| ≤ a − 1, t ∈ [0, T ], and using (6.3) we have Moreover, using (3.7) and (6.5) we have that for any γ = (t, ω) Let us take the supremum in γ ∈ T in the right-hand side and then fix T so small that LaT (L D Z 1 C( T * ) + 1) < 1/2.We obtain that Z 1 − Z 2 C( T ) ≤ Z 1 − Z 2 C( T ) /2, which leads to the contradiction Z 1 − Z 2 C( T ) = 0.The proof of the theorem is complete.
Proof of of Theorem 3.1 Now, we consider the specific system (2.3)-(2.7).The conditions (6.4)-(6.6)are apparently fulfilled in this case.Inspecting the proof of Theorem 6.1 we see that the conditions (6.2)-(6.3)are not used in parts 1-5 (they are only used in parts 6 and 7).We have proved (in parts 1-5) that for any D ∈ L ∞ (0, T ) with 0 ≤ D(t) ≤ a − 1 for a.e.t ∈ [0, T ], there exists a non-negative Lipschitz continuous solution Z [D] ∈ K T ,a of (3.7)-(3.8),hence of (2.4)-(2.7).The existence of Z [D] was obtained due to the contractivity of the operator F [D] on K T ,a .Hence, Z [D] can be considered as the uniform limit of the sequence of functions {Z k } generated by Notice that due to (6.4) and the definition of F [D] we have which implies the estimate . (6.15) Further, we shall choose T satisfying the last inequality, in addition to (6.7).

123
Denote by 0 the set points in (0, 1) on which S 0 or I 0 is non-differentiable, together with the points ω = 0 and ω = 1.Denote # := {γ ∈ T : ξ f (γ ) ∈ 0 or ξ g (γ ) ∈ 0 } = {ω f [(ξ, 0)](t), ω g [(ξ, 0)](t) : This set consists of finite number of curves in T .Observe that the assumption f , g < 0 on (0, 1) implies that the set ¯ T = T \ # consists of finite number of open sets, further called facets.We remind that γ f (γ ) and γ g (γ ) have Lipschitz derivatives in a neighborhood of T .Then the function Z 0 (γ ) = Z # (γ ) = ( S0 (γ f (γ )), Ī 0 (γ g (γ ))) is differentiable with a Lipschitz derivative on each facet of ¯ T \ # .In addition, for every γ ∈ T the functions s → σ (ω f [γ ](s)) and s → σ (ω g [γ ](s)) are differentiable and have Lipschitz derivatives on every of the finite number of intervals for s in which the argument of σ belongs to one facet.The same applies to the functions ρ and μ.Thanks to the properties mentioned in this paragraph, we can differentiate Z k+1 with respect to γ ∈ ¯ T using (6.10)-(6.11).Skipping the cumbersome details, we obtain the following relations: where Lip # (Q) is a common Lipschitz constant of a function Q on each facet of ¯ T (for Q : ¯ T → R 2 which is Lipschitz on every facet), Lip # dS 0 dω is the Lipschitz constant of dS 0 dω on each of the intervals of its existence, c 1 and c 2 are constants (which may depend on Lip(D) and Z 0 C ), c is a constant which is independent of Z 0 and D with 0 ≤ D ≤ a − 1.The derivation of this recurrent inequality also uses the fact that F [D] is an affine mapping of Z .Since Lip # dS 0 dω is finite, we obtain inductively that for every k (T c) j ≤ 2 c 1 + c 2 Lip # dS 0 dω , provided that 2cT ≤ 1.We add the last condition to (6.15) and (6.7).The choice of T is still independent of the initial distribution Z 0 and the particular D.

Fig. 1
Fig. 1 Characteristic lines and illustration of the notations

. 4 )
Vice versa, any pair of continuous functions (z S [•](•), z I [•](•)) defined on the sets as in the previous exposed lines determines a continuous pair (S, I ) : → R 2 by the relations

Fig. 2 Fig. 3
Fig. 2 Evolution of epidemiological population groups without vaccination

Fig. 4
Fig. 4 Evolution of epidemiological groups with no mortality rate

Fig. 6
Fig. 6 Saved lives per vaccine depending on the initial date of implementing vaccination and the vaccination rate v

Fig. 10
Fig. 10 The left plot presents efficient frontier for the three control scenarios; the red dots on the Pareto curves indicate the Pareto points corresponding to the baseline case.The right plot depicts the timedependence of the number of newly vaccinated individuals in the three control scenarios in the baseline case
7)-(3.8) has a Lipschitz continuous solution Z [D] = (S[D], I [D]) on T .Here T is independent of Z 0 .The Lipschitz constant of Z [D] can be estimated by the constant λ * := ( Z 0 The next step is to proof that the solution (S[D], I [D]) of (3.7)-(3.8)on [0, T ] depends in a Lipschitz way on D in a sense that will become clear in the next lines.For any two functions D For the limit Z [D] of Z k we obtain Lip # ∂ Z [D] ∂γ ≤ 2 c 1 + c 2 Lip # dS 0 dω .Since every horizontal and every vertical line intersects # only finite number of times, the partial derivatives of Z [D] in each of the variables (t, ω) exist, except a finite number of points, for every value of the other variable.the obtained differentiability properties of the solution (S[D], I [D]), we may employ D](t, ω) + I (t, ω)) dω = − D](t, ω) + I [D](t, ω)) dω ≥ 1 − μ L ∞ (0